% This file analyzes the experiments intended to reveal the tensor monopole
% generates Fig1CD
flag = {[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,-1,0],[1,0,1],[1,0,-1],[0,1,1],[0,1,-1]};

%% Fig1CD
load('./data/Fig1CD.mat')
alpha_fit=linspace(0,pi/2,41);
rabi_estimate_save = QGT_rabi_freq_detune(alpha_fit,0,0,0,2e6,flag)/1e6;
% Rabi oscillations
figure
% Fig1C

cc=[0.9647    0.3294    0.0510
    0.9020         0    0.0627
    0.0941    0.5294    0.9098
    0.8941    0.2314    0.5333
    0.1020    0.5686    0.4745
    0.9882    0.7098    0.0588
    0.5686    0.7490    0.1020
    0.0549    0.2706    0.6275
    0.0980    0.5412    0.1843];

% SQ
hold on
for hlp=1:length(flag)
    switch hlp
        case {4,6}
            plot(alpha_fit/pi*180,reshape(rabi_estimate_save(1,hlp,:),1,[]),'--','LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
        case {5,7}
            plot(alpha_fit/pi*180,reshape(rabi_estimate_save(1,hlp,:),1,[]),':','LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
        otherwise
            plot(alpha_fit/pi*180,reshape(rabi_estimate_save(1,hlp,:),1,[]),'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
    end
end
for hlp=1:length(flag)
    switch hlp
        case {1,2,3} % diagonal terms, solid circles
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),reshape(drabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
        %case {2,3} % diagonal terms, solid circles
        %    errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),reshape(drabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
        case {5}
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),reshape(drabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
        case {4,6}
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),reshape(drabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
        case {7,8}
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),reshape(drabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),'^','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
        case {9}
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),reshape(drabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),'^','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
        otherwise
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),reshape(drabi_fit_save(1,hlp,:),1,[])./reshape(E_pert_save(1,hlp,:),1,[]),'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
    end
end

% legend('\Gamma^\alpha','\Gamma^\beta','\Gamma^\phi','\Gamma^{\alpha\beta}','\Gamma^{\alpha\bar{\beta}}',...
%     '\Gamma^{\alpha\phi}','\Gamma^{\alpha\bar{\phi}}','\Gamma^{\beta\phi}','\Gamma^{\beta\bar{\phi}}','Interpreter','latex')

ymax=(max(max(rabi_estimate_save(1,:,:))))+0.2;
ylim([0,ymax]);xlim([0,90]);set(gca,'XTick',[0:30:90])
xlabel('\alpha(\circ)','FontSize',18)
ylabel('Matrix Element (MHz)','FontSize',18)
title('SQ','FontSize',20)
hold off
set(gca,'FontSize',18)

% subplot(1,2,2) % DQ
figure
hold on
for hlp=1:length(flag)
    switch hlp
        case {4,6}
            plot(alpha_fit/pi*180,reshape(rabi_estimate_save(2,hlp,:),1,[]),'--','LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
        case {5,7}
            plot(alpha_fit/pi*180,reshape(rabi_estimate_save(2,hlp,:),1,[]),':','LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
        otherwise
            plot(alpha_fit/pi*180,reshape(rabi_estimate_save(2,hlp,:),1,[]),'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
    end
end
for hlp=1:length(flag)
    switch hlp
        case {1,2,3} % diagonal terms, solid circles
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),reshape(drabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
        %case {2,3} % diagonal terms, solid circles
        %    errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),reshape(drabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
        case {5}
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),reshape(drabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
        case {4,6}
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),reshape(drabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
        case {7,8}
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),reshape(drabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),'^','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
        case {9}
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),reshape(drabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),'^','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
        otherwise
            errorbar(alpha_sweep/pi*180,reshape(rabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),reshape(drabi_fit_save(2,hlp,:),1,[])./reshape(E_pert_save(2,hlp,:),1,[]),'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
    end
end
ymax=(max(max(rabi_estimate_save(2,:,:))))+0.2;
ylim([0,ymax]);xlim([0,90]);set(gca,'XTick',[0:30:90])
title('DQ','FontSize',20)
hold off
set(gca,'FontSize',18)
xlabel('\alpha(\circ)','FontSize',18)
ylabel('Matrix Element (MHz)','FontSize',18)
lgd=legend('$$\Gamma_{\alpha}$$','$$\Gamma_{\beta}$$','$$\Gamma_{\phi}$$','$$\Gamma_{\alpha\beta}$$','$$\Gamma_{\alpha\overline{\beta}}$$',...
    '$$\Gamma_{\alpha\phi}$$','$$\Gamma_{\alpha\overline{\phi}}$$','$$\Gamma_{\beta\phi}$$','$$\Gamma_{\beta\overline{\phi}}$$');
set(lgd,'Interpreter','latex');lgd.FontSize=18;


function [rabi_g,E_total,full_map,g_munu_save] = QGT_rabi_freq_detune(alpha_sweep,beta,phi,rd,r,flag)

rabi_g = zeros(2,length(flag),length(alpha_sweep));
E_total = zeros(2,length(alpha_sweep)); % 1: SQ, 2: DQ
full_map = zeros(1,length(alpha_sweep));
dpert=0.005;
g_out = zeros(length(flag),length(alpha_sweep));
g_munu_save = zeros(6,length(alpha_sweep));

for hlp=1:length(alpha_sweep)
    theta=alpha_sweep(hlp);
    H0=Ham_detune(theta,beta,phi,rd);
    [V,E]=eig(H0);
    [E,idx]=sort(diag(E));
    
    E_total(:,hlp)=[abs(E(2)-E(1))*2;abs(E(3)-E(1))]*r; % note we already take 2X the SQ frequency here
    
    V=V(:,idx);
    um=V(:,1)/sqrt(V(:,1)'*V(:,1));
    u0=V(:,2)/sqrt(V(:,2)'*V(:,2));
    up=V(:,3)/sqrt(V(:,3)'*V(:,3));
    
    for hlp_flg=1:length(flag)
        flg_tmp=flag{hlp_flg};
        
        theta_oscillate=flg_tmp(1);
        beta_oscillate=flg_tmp(2);
        phi_oscillate=flg_tmp(3);
        
        for hlp_omega=1:2 % omega
            dH=Ham_detune(theta+dpert*theta_oscillate,beta+dpert*beta_oscillate,phi+dpert*phi_oscillate,rd);
            dH=(dH-H0)/dpert;
            
            switch hlp_omega
                case 1
                    rabi_g(hlp_omega,hlp_flg,hlp)=abs(um'*dH*u0*r);
                    
                    g_out(hlp_flg,hlp)=g_out(hlp_flg,hlp) + ...
                        abs(um'*dH*u0*r)^2/((E(1)-E(2))*r)^2;
                    
                case 2
                    rabi_g(hlp_omega,hlp_flg,hlp)=abs(um'*dH*up*r);
                    
                    g_out(hlp_flg,hlp)=g_out(hlp_flg,hlp) + ...
                        (rabi_g(hlp_omega,hlp_flg,hlp))^2/((E(1)-E(3))*r)^2;
            end
        end
    end
    
    g_munu_save(1,hlp)=g_out(1,hlp);
    g_munu_save(2,hlp)=g_out(2,hlp);
    g_munu_save(3,hlp)=g_out(3,hlp);
    g_munu_save(4,hlp)=(g_out(4,hlp)-g_out(5,hlp))/4;
    g_munu_save(5,hlp)=(g_out(6,hlp)-g_out(7,hlp))/4;
    g_munu_save(6,hlp)=(g_out(8,hlp)-g_out(9,hlp))/4;
    
end
end

function Ham_out = Ham_detune(alpha,beta,phi,rd) %
Ham_out=[0,cos(alpha)*exp(-1i*beta),0;cos(alpha)*exp(1i*beta),0,sin(alpha)*exp(1i*phi);...
    0,sin(alpha)*exp(-1i*phi),0]+diag([rd,0,-rd]);
end

function [QT,QT_err] = Herror_charge(alpha_sweep,H_berry_curvature, H_berry_err)
% calculates the topological charge
% takes alpha_sweep which is not equidistant

QT=sum((H_berry_curvature(1:end-1)+H_berry_curvature(2:end)).*(alpha_sweep(2:end)-alpha_sweep(1:end-1)));

QT_err=sqrt(sum((H_berry_err(1:end-1).^2+H_berry_err(2:end).^2).*(alpha_sweep(2:end)-alpha_sweep(1:end-1)).^2));
end

